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Abstract 

We revisit the question whether the critical behavior of sandpile models with 
sticky grains is in the directed percolation universality class. Our earlier theoretical 
arguments in favor, supported by evidence from numerical simulations [ Phys. Rev. 
Lett., 89 (2002) 104303], have been disputed by Bonachela et al. [Phys. Rev. E 74 
(2004) 050102] for sandpiles with no preferred direction. We discuss possible reasons 
for the discrepancy. Our new results of longer simulations of the one-dimensional 
undirected model fully support our earlier conclusions. 



After the pioneering work of Bak, Tang and Wiesenfeld[l] in 1987, sand- 
pile models have been studied extensively in the past two decades, both as 
paradigms of self-organized critical systems in general[2], and also as mod- 
els of real granular matter [3]. Many different types of sandpile models with 
different toppling rules have been studied [4] : deterministic and stochastic, 
with or without preferred direction, different instability criteria [5], or particle 
distribution rules [6], with fixed energy [7] etc.. Most of these models could 
only be studied numerically, and for a while it seemed that each new variation 
studied belonged to a new universality class of critical behavior. Though not 
complete, a broad picture of the different universality classes of self-organized 
critical behavior has emerged in recent years [8,9]. 

In an earlier paper [10], we have argued that the generic behavior of sandpile 
models is in the universality class of directed percolation (DP), and mod- 
els with deterministic toppling rules like the original BTW model, and mod- 
els with stochastic toppling rules like the Manna models, are unstable to a 
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perturbation of introduction of "stickiness" in the toppling rules, and under 
renormalization. the flows are directed towards the DP-fixed point. These ar- 
guments are reasonable, but not rigorous, and we presented a detailed study 
of a specific model, where some of the steps in the arguments could be shown 
to be vahd, and we used detailed Monte Carlo simulations to check our con- 
clusions. 

Some of the conclusions of this paper have recently been disputed by Bonachela 
et al. [11]. These authors contend that while the directed sandpiles with sticky 
grains show DP behavior in the SOC limit, our arguments do not apply to 
undirected sandpiles, where, even with stickiness, the critical behavior contin- 
ues be the same as that of the Manna model, i. e., in the universality class of 
the directed percolation with a conservation law (hereafter referred to as the 
Manna/C-DP universality class). 

In this paper, we will try discussing these conflicting claims, and also present 
some data from more recent extensive simulations, which supports our original 
conclusions. We start by defining the model precisely, and then summarize the 
arguments of [10]. We then discuss the simulations of [11], and finally present 
the results of our new more extensive simulations. 

First the precise definition of the model. We consider the directed model on a 
(1 + l)-dimensional square lattice, for definiteness. Generalizations to higher 
dimensions are straight forward. The sites on an L x M torus are labelled by 
euclidean coordinates {i-ij) with (^ + j) even and j increasing downward. At 
each site there is a non-negative integer hij to be called the height of 

the pile at that site. Initially all hi^j are zero. The system is driven by choosing 
a site at random and increasing the height at that site by one. 

The 'stickiness' of the grains is characterized by a parameter p, and its role 
in the dynamics of sandpiles is defined as follows: A site is said to become 
unstable at time t, if at least one particle is added to it at time t, and its 
height becomes greater than 1. A site made unstable at time t relaxes 

at the time {t + 1) stochastically: With probability (1 — p), it becomes stable 
without losing any grains, and the added particle (s) sticks to the existing 
grains. Otherwise (with probability p), the site topples, and the height at 
the site decreases by two, and the site becomes stable. We introduce bulk 
dissipation: at each toppling, with probability 6 both grains from the toppling 
are lost, other wise (with probability 1 — 5), the two grains are transferred to 
the two downward neighbors {i ± 1, j + 1). 

We relax all the unstable sites by parallel dynamics. An unstable site is relaxed 
in one step, independent of whether it received one or more grains at the 
previous time step. Once a site has relaxed, it remains stable until perturbed 
again by new grains coming to the site. This relaxation process is repeated 
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until all sites become stable, and then a new grain is added. 

The model is specified by two parameters p and 5. For this model, the follow- 
ing results can be proved [10]: 

i) Depending on the values of p and 6, two different behaviors are possible. 
There is a threshold p*{6), such that for p > p*{d), there is a steady state 
in the system, but for p < p* {6) , no steady state is possible, and the mean 
height of the pile grows linearly with time [Fig. 1]. 

ii) The boundary linep = p*{6) is exactly given in terms of the function n_Dp(A), 
which gives the mean number of infected and boundary sites in a cluster in 
a directed site percolation process on the same lattice with infection prob- 
ability A. In particular, this boundary line meets the S — line at p — Pc, 
where Pc is the exact directed site percolation threshold for the square lat- 
tice. 

iii) At the boundary line p = p*{S), the distribution of sizes of avalanches is ex- 
actly the same as in the DP process, with a infection probability A = p{l—5). 

iv) In the regime where the steady state exists, mean cluster size is finite for 
non-zero S, and we get SOC behavior only ior pc < p < 1, and 5 — > 0. 

v) The model can be solved exactly along three of the boundary fines, p — 1, 
S — 1 and p — 0. 

We find that correlations between heights at different sites in the steady state 
are quite weak, and the steady state is almost a product measure state. We ar- 
gued that if these correlations are irrelevant, one can study the avalanche pro- 
cess in a background where these correlations arc absent. Then the avalanche 
process becomes a Domanay- Kinzel process with two parameters {pi,P2), 
where P2 = p and pi is determined in terms of the concentration of sites with 
height zero in the steady state. If the correlations in height are indeed irrele- 
vant, the avalanche distribution function for all p along the SOC line in Fig. 
1 would show DP exponents, except the end point p — 1,5 — where the 
model has deterministic toppling rules. 

We then extended the arguments given to undirected sandpiles, by considering 
the time-evolution of a d-dimensional pile as a infection process on a ((i -|- 1)- 
dimensional lattice. It can be proved that point i)-iv) listed above continue 
to hold for the undirected case, and p*{5 = 0) is exactly given by the critical 
threshold of the corresponding DP process. 

Again, if we can neglect correlations between heights of pile at different sites 
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Fig. 1. The phase diagram in the p-6 plane. The inset shows the variation of 
[pc — P*{S)] versus S, with the straight hne showing the theoretical fit (from [10]). 

in the steady state, the avalanche exponents in the model can be seen to be 
same as DP exponents. Correlations between heights in this case are stronger 
than in the directed height at a site remains unchanged if there is no 

activity, and is obviously correlated in time. Thus, it is not clear that these 
are also irrelevant in the renormalization group sense in our problem. In fact, 
Bonachcla ct al argue that these arc not. However, we note that because our 
toppling rules do not depend on the height, so long as it is greater than 1, 
most of these correlations do not come into picture, as probabilities of different 
topplings depend only whether the height at a site was zero or non-zero before 
the addition of the new particle that made it unstable. The agreement of the 
results of our numerical simulations with the theoretical predictions shows 
that indeed, for determining the critical exponents, neglecting the effect of 
these correlations is justified. 

Bonachela et al. have argued that while the arguments in [10] are correct for 
the directed model, the neglect of correlations is not valid in the case of undi- 
rected models. They presented Monte Carlo evidence based on simulations of 
the fixed-energy sandpile (FES) model, and also some nonrigorous arguments 
in support of their proposition. 

Let us discuss the simulations of the FES first. It should be noted that one 
cannot get any avalanche exponents directly from such simulations, as by 
definition, the FES shows a single avalanche that never stops, and there is 
no distribution of avalanche sizes to quantify. The avalanche exponents can 
only be inferred indirectly, by determining some other exponents, and then 
using scaling theory to relate them to the conventional avalanche exponents. 
For example, exponent governing the decay of activity-activity correlation 
function with time in the steady state of the fixed-energy sandpile is used to 
determine the exponent 
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Fig. 2. Schematic representation of time-evolution in the model in different cases: 
(a) non-overlapping avalanche clusters for very small e (b) an infinite avalanche 
cluster and disconnected finite clusters for intermediate e, (c) Avalanche cluster in 
the fixed energy sandpile. 

The FES sandpile model can be obtained from our model by considering a 
uniform addition of particles at rate e per site per unit time, with e <^ 5, and 
take the double limit e — 0, and 6 ^ 0. In the usual slowly driven sandpiles, 
one takes the limit e first, and then 5 — limit. In the FES case, we 
have to take the limit e — > with mean activity d = e/{26) held fixed, and 
then take the limit a 0. The different order of limits can lead to different 
scaling behaviors. 

In Fig. 2, we have shown a schematic space-time plot of the evolution of the 
activity in a part of the sandpile in different cases (The y-axis is time, and it 
increases upwards). In the case of slowly driven sandpile case , with e ^ 5, we 
see that different activity clusters are disconnected from each other, and show a 
wide variation of sizes [ panel (a)]. If e is comparable to 6, there are overlapping 
activity clusters [panel (b)]. However, there still are disconnected clusters of 
activity, and one can find events where activity starts from inert background ( 
due to addition of a particle), and the avalanche activity so generated diffuses 
around, and may later die, or merge with the infinite avalanche. In (c), we 
show the FES activity cluster. There are no disconnected activity clusters, 
and the large-scale structure of the cluster is different from (b), in that no 
activity can start from inert state unless a neighbor was active. Thus, the 
large-scale structure of clusters is clearly different in the three cases. 

In particular, we note that the avalanche distributions in the conventional 
slowly driven sandpile are determined by the properties of disconnected finite 
clusters (when e -C (5, the probability that the added new particle will be added 
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to an active site is small). In contrast, the FES properties are determined by 
the properties of the single infinite cluster. 

We note, as argued in [10], that on the fine p = p*{6), the probabilities of all 
clusters in the sticky sandpile model are exactly those of DP clusters, with 
the infection probability A given as an explicit function of p and S. This is so, 
as on this line, the fraction of sites with height is zero in the steady state, 
and every infected site has the same probability A = p(l — 6) of infecting 
a downward neighbor. If p > p*{6) by a small amount Ap, the mean height 
would be large, and concentration of sites with height in the steady state is 
small of order Ap. The size of a cluster s is defined in the space time plot as the 
total number of sites whose height has been altered during a single avalanche. 
Clearly the difference in cluster size distribution |Prob(s) — Prob£)p(s)| tend 
to zero as Ap tends to zero, and the probabihties for small avalanches tend to 
the DP value near the boundary p = P*{S). 

We also note the activity cluster in the FES for large average energy E is the 
same as the infinite DP-cluster with site percolation parameter p, as each site 
having at least one particle coming to it would topple with a probability p. 
Thus, in the limit of large mean energy, all n-point correlation functions of 
the activity in the steady state are exactly the same as in the DP process with 
infection parameter A = p. 

The only possibility for seeing the Manna exponents in the avalanche statistics 
then is to look for it in some intermediate range of values of s , Smin ^ s ^ 
Smax, where Smin is the crossover scale above which DP behavior is presumably 
lost, and Smax is the upper cutoff imposed by finite value of 6. Here Smin is 
expected to diverge as {Ap}^"-' near the critical line. As the difference between 
the DP and C-DP exponents is small, Smin would be expected to be large, 
and the difference (even if present) is difficult to establish convincingly by 
simulations. 

Bonachcla et al. have also studied the problem using a numerical study of the 
coupled Langevin equations for the density fields for activity and particles. 
However, this study is also a numerical integration of the stochastic evolution 
equations, and is qualitatively not different from a Monte Carlo simulation 
working with coarse-grained variables instead of the original variables defined 
on the lattice sites. 

Finally, we discuss our more recent Monte Carlo data for the test case of 
one-dimensional undirected sandpile model with sticky grains. We monitored 

the average transverse cluster size as a function of the number of topplings 
in the cluster. This is expected to be less sensitive to the upper cutoff on 
the s. We did our simulations for p = 0.85, 6 = 10~^ on a line of length 
L — 4096 with periodic boundary conditions. This value of 5 is a factor 10 
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Fig. 3. Log-log plot of mean transverse size of cluster as a function of number 
of topplings in the cluster for the sandpile model with sticky grains(SG) in one 
dimension for p = 0.85 (blue line). For a comparison, we also plot simulation results 
for DP clusters (red line), and Manna model clusters (green line). 

smaller than the values used in [10], which allows us to generate much larger 
clusters. The chosen value of p (and a small 5) is half-way between and 
1, and one would expect the effects of crossovers from DP and deterministic 
limits to be small. The mean height in the steady state at this value of p is 
1.977(5), and approximately 8.3% of the sites are of zero height. Thus, any 
deviations from the DP behavior coming from the presence of such sites should 
be measurable. The results of data taken for 2 x 10^ avalanches are shown in 
Fig. 3. For comparison, we have also plotted the results for DP clusters with 
p = 0.7, and for the Manna model. It is clear that the observed slope is much 
closer to DP than to that for Manna clusters. 

The theoretical arguments for neglecting, or not neglecting the correlations are 
not fully convincing, and appeal to numerical simulations for their justifica- 
tion. Extracting critical exponents from simulation data could be complicated 
by crossover effects, as the difference in the avalanche exponents for the Manna 
and DP cases is not large. We hope that further work on this problem will 
clarify the situation, and lead to better insight into the problem. 
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